楕円描画について

%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib.pyplot import *

import numpy as np

pi=np.pi
sqrt=np.sqrt
cos = np.cos
sin = np.sin

1.傾いた楕円の描画方法

誤差楕円の描画をしたいときがあるのでメモ

普通の楕円

基本の楕円 $$\frac{1}{4} x_1^2 + x_2^2 =1$$

描画は媒介変数表示が便利． $$\begin{eqnarray} x_1 &=& 2 \cos \theta \\ x_2 &=& \sin \theta \end{eqnarray}$$

上記の式を$\theta=[-\pi,\pi]$の範囲内で適当なステップで区切ればok

一般化すると $$ax_1^2 + bx_2^2 =1$$ に対して媒介変数表示は $$\begin{eqnarray} x_1 &=& \frac{1}{\sqrt{a}} \cos \theta \\ x_2 &=& \frac{1}{\sqrt{b}} \sin \theta \end{eqnarray}$$

全体的に拡大・縮小するなら半径を変える感じで $$ax_1^2 + bx_2^2 =\chi^2$$ に対して媒介変数表示は $$\begin{eqnarray} x_1 &=& \frac{\chi}{\sqrt{a}} \cos \theta \\ x_2 &=& \frac{\chi}{\sqrt{b}} \sin \theta \end{eqnarray}$$

""" parameter """    
a,b,chi2=1/4,1,2

def abc2elip(a,b,chi2=1,resol=100): #ax+by=chi^2
    th=np.linspace(-pi,pi,resol)
    x=np.c_[ sqrt(chi2/a)*cos(th),sqrt(chi2/b)*sin(th)].T
    return x

x=abc2elip(a,b,chi2=chi2)
figure()
plot(x[0,:],x[1,:])
grid()

楕円を傾ける

おなじみの回転行列で傾ける． 角度$\theta$だけ回転させた後の点を$\boldsymbol{x}'=[x_1',x_2']^\mathrm{T}$とおくと，以下の変数変換でok $$\left( \begin{array}{c} x_1' \\x_2' \end{array} \right) =\left( \begin{array}{c} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array} \right) \left( \begin{array}{c} x_1 \\x_2 \end{array} \right)$$

""" parameter """    
th=pi/6

def abcth2elip(a,b,th,chi2=1,resol=100):# th=radian
    x=abc2elip(a,b,chi2,resol)
    R=np.array([[cos(th),-sin(th)],[sin(th),cos(th)]])
    return R@x

x=abcth2elip(a,b,th,chi2=chi2)
figure()
plot(x[0,:],x[1,:])
grid()

二次形式から描画する

基本楕円

基本の楕円は以下の式で表すことができる．(二次形式という) $$\begin{eqnarray} \left( \begin{array}{c} x1 \\x2 \end{array} \right) ^\mathrm{T} \left( \begin{array}{c} \frac{1}{4}& 0 \\ 0&1 \end{array} \right) \left( \begin{array}{c} x1 \\x2 \end{array} \right) &=&1\\ \equiv \boldsymbol{x}^\mathrm{T}\boldsymbol{Ax}&=&1 \end{eqnarray}$$

楕円の表現行列$\boldsymbol{A}$が対角行列のときは以下の媒介変数表示でok $$\begin{eqnarray} x_1 &=& \frac{1}{\sqrt{A(1,1)}} \cos \theta \\ x_2 &=& \frac{1}{\sqrt{A(2,2)}} \sin \theta \end{eqnarray}$$

""" parameter """    
Vinv=np.array([[a,0],[0,b]])
print(Vinv)

[[0.25 0.  ]
 [0.   1.  ]]

def ivar2elip(Vinv,chi2=1,resol=100): # (x^T)Qx=chi2  V=diag(var_x1,var_x2) //Quadratic form
    x=abc2elip(Vinv[0,0],Vinv[1,1],chi2,resol)
    return x

x=ivar2elip(Vinv,chi2=chi2)
figure()
plot(x[0,:],x[1,:])
grid()

傾いた楕円

表現行列$\boldsymbol{A}$が対角行列でないとき，傾いた楕円になる． このパターンが一番多いだろうし厄介．

まず表現行列$\boldsymbol{A}$の固有値問題を解く． 解いた結果，固有値$\lambda_1,\lambda_2$，固有ベクトル$\boldsymbol{p}_1,\boldsymbol{p}_2$が得られる． 固有ベクトルを横に並べてできる行列を$\boldsymbol{P}$と置く．このとき，楕円の媒介変数表示は以下で得られる．

$$\left( \begin{array}{c} x_1 \\ x_2 \end{array} \right)= \boldsymbol{P} \left( \begin{array}{c} y_1 \\ y_2 \end{array} \right)$$

ここで， $$\begin{eqnarray} y_1 &=& \frac{1}{\sqrt{\lambda_1}} \cos \theta \\ y_2 &=& \frac{1}{\sqrt{\lambda_2}} \sin \theta \end{eqnarray}$$

""" parameter """    
def compose(Vinv,th): # get representation matrix
    R=np.array([[cos(th),-sin(th)],[sin(th),cos(th)]])
    return R@Vinv@R.T

def decompose(Qinv): # representation mat ->  Vinv,Rot(th)
    lam,R=np.linalg.eig(Qinv)
    Vinv=np.diag(lam)
    return (Vinv,R)

Qinv=compose(Vinv,th)
print(Qinv)

[[ 0.4375     -0.32475953]
 [-0.32475953  0.8125    ]]

def icov2elip(Qinv,chi2=1,resol=100): # (x^T)Qx=chi2  Q=Cov(x)//Quadratic form
    Vinv,R=decompose(Qinv)
    x=ivar2elip(Vinv,chi2,resol)
    return R@x

x=icov2elip(Qinv,chi2=chi2)
figure()
plot(x[0,:],x[1,:])
grid()

理論武装編

対角行列$\varLambda$を表現行列とするような基本の楕円を考える． $$\boldsymbol{x}^\mathrm{T} \boldsymbol{\varLambda x}=1$$

適当な回転行列$\boldsymbol{R}$によって基本の楕円を回転させる． $$\boldsymbol{x}'=\boldsymbol{Rx}$$

上記の変数変換によって二次形式表現は以下のようになる． $$\begin{eqnarray} (\boldsymbol{R}^{\mathrm{T}}\boldsymbol{x}')^\mathrm{T} \boldsymbol{\varLambda} (\boldsymbol{R}^{\mathrm{T}}\boldsymbol{x}')&=&1 \\ \iff \boldsymbol{x}'^{\mathrm{T}} \boldsymbol{R} \boldsymbol{\varLambda} \boldsymbol{R}^{\mathrm{T}} \boldsymbol{x}' &=&1 \\ \equiv \boldsymbol{x}'^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}' &=&1 \end{eqnarray}$$

ここで，行列$\boldsymbol{A}$は傾いた楕円の表現行列である．表現行列の中身を見ると，(当たり前だが)基本の楕円の表現行列と，回転行列という2つの要素から構成されている． ということは，逆に，傾いた楕円の表現行列を，基本楕円の表現行列と回転行列に分解可能なはずだ，というモチベーションが生まれる．

そこで行列$\boldsymbol{A}$を対角化してみる． 対角化には，まず行列$\boldsymbol{A}$の固有値問題を解く． そして，固有値$\lambda_1,\lambda_2$,固有ベクトル$\boldsymbol{p}_1,\boldsymbol{p}_2$を用いて，以下の対角行列$\boldsymbol{\Sigma}$，直交行列$\boldsymbol{P}$を作る． $$\begin{eqnarray} \boldsymbol{\Sigma}&=& \left( \begin{array}{c} \lambda_1& 0 \\ 0&\lambda_2 \end{array} \right) \\ \boldsymbol{P}&=& \left(\boldsymbol{p}_1,\boldsymbol{p}_2 \right) \end{eqnarray}$$

最後に上記の行列を用いて以下のように対角化が行われる． $$\boldsymbol{P}^{\mathrm{T}} \boldsymbol{AP}=\boldsymbol{\Sigma}$$

$\boldsymbol{P}^{\mathrm{T}} \boldsymbol{AP}=\boldsymbol{\Sigma}$と$\boldsymbol{R \varLambda} \boldsymbol{R}^{\mathrm{T}}=\boldsymbol{A}$を比較する． ここは論理が飛躍するが，ぱっと見で$\boldsymbol{\Sigma}=\boldsymbol{\varLambda}$，$\boldsymbol{P}=\boldsymbol{R}$が成り立つと思われる．

(恐らく，$\boldsymbol{R} \boldsymbol{\varLambda} \boldsymbol{R}^{\mathrm{T}} = \boldsymbol{A}$の関係式より，$\boldsymbol{A}$と$\boldsymbol{\varLambda}$が相似である(多分)ので，相似な行列の固有値は等しいことを利用して，$\boldsymbol{\Sigma}=\boldsymbol{\varLambda}$を導くことができる)

まとめると，傾いた楕円の表現行列について以下が成り立つ． $$\boldsymbol{P}^{\mathrm{T}} \boldsymbol{AP}=\boldsymbol{\varLambda}$$

したがって，$\boldsymbol{x}'=\boldsymbol{Px}$なる変数変換によって，傾いた楕円の二次形式表現は以下のように，基本楕円へと変換可能である． $$\begin{eqnarray} \boldsymbol{x}'^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}' &=&1 \\ \iff (\boldsymbol{Px})^{\mathrm{T}} \boldsymbol{A} (\boldsymbol{Px}) &=&1 \\ \iff \boldsymbol{x}^{\mathrm{T}} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{Px} &=&1 \\ \iff \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\varLambda} \boldsymbol{x} &=&1 \end{eqnarray}$$

ここで，基本楕円の表現行列の対角成分は，$\boldsymbol{\Sigma}=\boldsymbol{\varLambda}$より固有値なので，以下のようにパラメータ表示できる． $$\begin{eqnarray} x_1 &=& \frac{1}{\sqrt{\lambda_1}} \cos \theta \\ x_2 &=& \frac{1}{\sqrt{\lambda_2}} \sin \theta \end{eqnarray}$$

最後に，$\boldsymbol{x}'=\boldsymbol{Px}$なる変数変換を用いてパラメータ表示した基本楕円を傾ければ描画完了．

Vinv=np.array([[a,0],[0,b]])
print(Vinv)
Qinv=compose(Vinv,th)
print(Qinv)
Vinv,R=decompose(Qinv)
print(Vinv)
print(R)

[[0.25 0.  ]
 [0.   1.  ]]
[[ 0.4375     -0.32475953]
 [-0.32475953  0.8125    ]]
[[0.25 0.  ]
 [0.   1.  ]]
[[-0.8660254  0.5      ]
 [-0.5       -0.8660254]]

追記

描画したいだけならcontourを使って陰関数の描画をしたほうが楽

http://stackoverflow.com/questions/2484527/is-it-possible-to-plot-implicit-equations-using-matplotlib

http://stackoverflow.com/questions/4680525/plotting-implicit-equations-in-3d

def icov2implicit2(Q,chi2=1):
    return lambda x1,x2 : Q[0,0]*(x1**2)+Q[1,1]*(x2**2) \
                     +2*Q[0,1]*(x1*x2)+ -chi2

def icov2implicit3(Q,chi2=1):
    return lambda x1,x2,x3 : Q[0,0]*(x1**2)+Q[1,1]*(x2**2)+Q[2,2]*(x3**2) \
                     +2*Q[0,1]*(x1*x2)+2*Q[1,2]*(x2*x3)+2*Q[2,0]*(x3*x1) -chi2

def plot_implicit(fn,Chi=0,xrange=[-1,1],yrange=[-1,1],resol=100,label='__nolabel'):
    """
    plot implicit function fn(x,y)=Chi (Chi=const.)
    """
    x=np.linspace(xrange[0],xrange[1],resol)
    y=np.linspace(yrange[0],yrange[1],resol)
    x,y =np.meshgrid(x,y)
    z=fn(x,y)
    if label =='__nolabel':
        plt.contour(x,y,z-Chi,[0])
    else:
        plt.contour(x,y,z-Chi,[0],label=label)
        
figure()    
f=icov2implicit2(Qinv)
plot_implicit(f,xrange=[-3,3],yrange=[-3,3])
grid()

2. 多次元楕円体

多次元楕円体は3次元までなら立体として描画できるが，それ以上の次元になると立体描画できない． そこで，多次元楕円体を任意の平面に正射影することを考える．正射影すれば，n次元楕円体から $_nC_2$個の二次元グラフが得られ，楕円体を想像することができる．

本ノートでは例として3次元楕円体をx-y,y-z,x-z平面へそれぞれ射影したグラフを描くが，(恐らく)4次元以上の楕円体にも適用可能である．

"""elipsoid define"""
def compose3d(Vinv,thx,thy,thz):
    Rx=np.array([[1,0,0],[0,cos(thx),-sin(thx)],[0,sin(thx),cos(thx)]])
    Ry=np.array([[cos(thy),0,-sin(thy)],[0,1,0],[sin(thy),0,cos(thy)]])
    Rz=np.array([[cos(thz),-sin(thz),0],[sin(thz),cos(thz),0],[0,0,1]])
    R=Rz@Ry@Rx
    return R@Vinv@R.T

from mpl_toolkits.mplot3d import axes3d
a,b,c=1/4,1,3
chi2=2
th,phi,psi=pi/6,pi/6,pi/6
Vinv3=np.array([[a,0,0],[0,b,0],[0,0,c]])
Qinv3=compose3d(Vinv3,th,phi,psi)
# implicit form
fv=icov2implicit3(Vinv3,chi2=chi2)
fc=icov2implicit3(Qinv3,chi2=chi2)

基本楕円体

基本楕円体の$x_\mathrm{i} , x_\mathrm{j}$平面への射影は，以下の表現行列に対する二次元楕円の描画を行えば良い．

$$\boldsymbol{A}_{\mathrm{i,j}} = \left( \begin{array}{c} \boldsymbol{A}[i,i] &\boldsymbol{A}[i,j] \\\boldsymbol{A}[j,i] & \boldsymbol{A}[j,j] \end{array} \right)$$

def ivar2elipn(Vinv,chi2=1,resol=100): # return prj(i,j):2D projection to (x_i,x_j) plane
    return lambda i,j : abc2elip(Vinv[i,i],Vinv[j,j],chi2,resol)

##1 2d check
prj=ivar2elipn(Vinv,chi2=chi2)
x = prj(0,1)
figure()
plot(x[0,:],x[1,:])
grid()

def plot_implicit3d(fn, ax,bbox=(-1,1)):
    ''' create a plot of an implicit function
    fn  ...implicit function (plot where fn==0)
    bbox ..the x,y,and z limits of plotted interval
    reference : http://stackoverflow.com/questions/4680525/plotting-implicit-equations-in-3d    
    '''
    xmin, xmax, ymin, ymax, zmin, zmax = bbox*3
#    fig = plt.figure()
#    ax = fig.add_subplot(111, projection='3d')
    A = np.linspace(xmin, xmax, 100) # resolution of the contour
    B = np.linspace(xmin, xmax, 15) # number of slices
    A1,A2 = np.meshgrid(A,A) # grid on which the contour is plotted

    for z in B: # plot contours in the XY plane
        X,Y = A1,A2
        Z = fn(X,Y,z)
        cset = ax.contour(X, Y, Z+z, [z], zdir='z')
        # [z] defines the only level to plot for this contour for this value of z

    for y in B: # plot contours in the XZ plane
        X,Z = A1,A2
        Y = fn(X,y,Z)
        cset = ax.contour(X, Y+y, Z, [y], zdir='y')

    for x in B: # plot contours in the YZ plane
        Y,Z = A1,A2
        X = fn(x,Y,Z)
        cset = ax.contour(X+x, Y, Z, [x], zdir='x')
        
    
##2 3d check
prj=ivar2elipn(Vinv3,chi2=chi2)
x12,x13,x23 =prj(0,1),prj(0,2),prj(1,2)

p=5
fig = plt.figure()
ax = fig.gca(projection='3d')       
plot_implicit3d(fv,ax,bbox=(-p,p))
ax.set_zlim3d(-p,p)
ax.set_xlim3d(-p,p)
ax.set_ylim3d(-p,p)
ax.set_xlabel('x')
ax.set_xlabel('y')
ax.set_xlabel('z')

ax.plot(x12[0,:],x12[1,:],'m',zs=-3,zdir='z',lw=5)
ax.plot(x23[0,:],x23[1,:],'c',zs=-3,zdir='x',lw=5)
ax.plot(x13[0,:],x13[1,:],'y',zs=-3,zdir='y',lw=5)

C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:17: UserWarning: No contour levels were found within the data range.
C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:23: UserWarning: No contour levels were found within the data range.
C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:28: UserWarning: No contour levels were found within the data range.

[<mpl_toolkits.mplot3d.art3d.Line3D at 0x22ab29e1408>]

一般楕円体

三次元楕円体の表現行列を$\boldsymbol{A}$とする．

まず，この三次元楕円体をどの平面に射影するか？という場合分けをなくすために，ベクトル・行列の並び替えを行う．

いま，三次元楕円体を$x_1,x_3$平面に射影したいとする．このとき，三次元楕円体は以下のように表現する．

$$\left( \begin{array}{c} x_1 \\ x_3\\x_2 \end{array} \right)^\mathrm{T} \left( \begin{array}{c} \boldsymbol{A}[1,1] & \boldsymbol{A}[1,3]& \boldsymbol{A}[1,2]\\ \boldsymbol{A}[3,1] & \boldsymbol{A}[3,3]&\boldsymbol{A}[3,2]\\ \boldsymbol{A}[2,1]&\boldsymbol{A}[2,3]&\boldsymbol{A}[2,2] \end{array} \right) \left( \begin{array}{c} x_1 \\ x_3\\x_2 \end{array} \right)$$

上記の表現行列を$\boldsymbol{H}$,ベクトルを$\boldsymbol{q}$とおく． このような並べ替えにより，表現行列$\boldsymbol{H}$で表される楕円体を$q_1,q_2$平面に射影する問題に帰着できる．

def swapr(T,i,j):
        sq=np.arange(0,T.shape[0])
        idi,idj=sq==i,sq==j
        sq[idi],sq[idj]=j,i
        return T[sq,:]    
    
def swapc(T,i,j):
    return swapr(T.T,i,j).T

def swap(T,i,j):
    return swapc(swapr(T,i,j),i,j)

T=np.array([[11,12,13],[21,22,23],[31,32,33]])
print(T)
print(swap(T,1,2)) # swap T[i,:]<->T[j,:],T[:,i]<->T[:,j] from basic.py

def perm2d(T,n,m):# Covariance Matrix swapping 
    '''
    T[n,n]->T[0,0],T[m,m]->T[1,1] by sequencial swapping
    
    A=[11 12 13 ; 21 22 23 ; 31 32 33] 
    -> perm2d(A,1,2)=[11 13 12 ; 31 33 32 ; 21 23 22]
    '''
    T2=swap(T,0,0)#copy    
    n,m=min(m,n),max(n,m)
    for it in np.arange(n,0,-1):
        T2=swap(T2,it-1,it)
    for it in np.arange(m,1,-1):
        T2=swap(T2,it-1,it)
    return T2

print(perm2d(T,0,2))

[[11 12 13]
 [21 22 23]
 [31 32 33]]
[[11 13 12]
 [31 33 32]
 [21 23 22]]
[[11 13 12]
 [31 33 32]
 [21 23 22]]

表現行列$\boldsymbol{H}$で表される楕円体を$q_1,q_2$平面に射影して得られる楕円の表現行列$\boldsymbol{H}_\mathrm{pj}$は以下の式で与えられる．

$$\boldsymbol{H}_\mathrm{pj} = \boldsymbol{A} - \boldsymbol{B} \boldsymbol{D}^{-1}\boldsymbol{C}$$ $$\boldsymbol{H}=\left( \begin{array}{c} \boldsymbol{A}&\boldsymbol{B}\\\boldsymbol{C}&\boldsymbol{D} \end{array} \right)$$

ただし，$\boldsymbol{A}\in \mathbb{R}^{2\mathrm{x}2}$

def mat_prj(Qinv,i,j):# representation mat of projected to i,j plane
    Q=perm2d(Qinv,i,j)# Covariance Matrix swapping
    A,B,C,D=block4(Q,1,1) # Q -> [A B ; C D] A=2x2 matrix
    Qpt = A-B@ np.linalg.solve(D,C) # 
    return Qpt
    
def block4(T,i,j): #
    '''
    T -> [A B ; C D] A.shape=(i,j)
    '''
    A=T[0:i+1,0:j+1]
    B=T[0:i+1,j+1:]
    C=T[i+1:,0:j+1]
    D=T[i+1:,j+1:]
    return (A,B,C,D)

def icov2elipn(Qinv,chi2=1,resol=100): # return prj(i,j):2D projection to (x_i,x_j) plane
    def prj(i,j):
        Qp=mat_prj(Qinv,i,j) # representation mat of projected to i,j plane
        return icov2elip(Qp,chi2=chi2,resol=resol)    
    return prj

##1 2d check
prj=icov2elipn(Qinv,chi2=2)
x = prj(0,1)
figure()
plot(x[0,:],x[1,:])
grid()
chi2=2

prj=icov2elipn(Qinv3,chi2=chi2)
x12,x13,x23=prj(0,1),prj(0,2),prj(1,2)

fig = plt.figure()
ax = fig.gca(projection='3d')       
plot_implicit3d(fc,ax,bbox=(-p,p))#
ax.plot(x12[0,:],x12[1,:],'m',zs=-3,zdir='z',lw=5)
ax.plot(x23[0,:],x23[1,:],'c',zs=-3,zdir='x',lw=5)
ax.plot(x13[0,:],x13[1,:],'y',zs=-3,zdir='y',lw=5)
# ax.set_aspect('equal')
ax.set_zlim3d(-p,p)
ax.set_xlim3d(-p,p)
ax.set_ylim3d(-p,p)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:17: UserWarning: No contour levels were found within the data range.
C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:23: UserWarning: No contour levels were found within the data range.
C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:28: UserWarning: No contour levels were found within the data range.

Text(0.5, 0, 'z')

理論武装編

表現行列$\boldsymbol{H}$で表される楕円体を$q_1,q_2$平面に射影する問題をもう少し具体的にする．

$q_1,q_2$平面に射影された楕円とは，楕円体の法線ベクトル$\boldsymbol{n}$が，$\boldsymbol{e}_3,\boldsymbol{e}_4,\cdots$に直交するような，$q_1,q_2$の集合である．ただし，$\boldsymbol{e}_i$は単位ベクトルである．

楕円体の法線ベクトルは$\boldsymbol{Hq}$で与えられる．したがって，上記の直交条件は以下のように表される．

$$\begin{eqnarray} \boldsymbol{e}_3^\mathrm{T} \boldsymbol{Hq} &=&0 \\ \boldsymbol{e}_4^\mathrm{T} \boldsymbol{Hq} &=&0 \\ \vdots \\ \iff \boldsymbol{E}^\mathrm{T}\boldsymbol{Hq} &=& \boldsymbol{0} \end{eqnarray}$$

ここで，$\boldsymbol{E}=[\boldsymbol{e}_3,\boldsymbol{e}_4,\cdots]$である．

三次元の場合は以下が得られる． $$(\boldsymbol{H}[1,1] q_1+\boldsymbol{H}[2,2]q_2)+\boldsymbol{H}[3,3]q_3=0$$

直交条件を使って楕円体方程式から$q_1,q_2$以外の要素を消去したい．

ということは，以下のことをしなくてはいけない．

1. 直交条件を$f(q_3,q_4,\cdots) = g(q_1,q_2)$の形にする．
2. 楕円体方程式をどうにかして$h(q_1,q_2) + f(q_3,q_4,\cdots) = const.$に変形する．
3. 楕円方程式に直交条件を代入し，$h(q_1,q_2) + g(q_1,q_2) = const.$を得る．

1.の手順について,3次元の場合からn次元を類推してみる．

三次元の場合は以下のように直交条件が変形できる． $$\begin{eqnarray} (\boldsymbol{H}[1,1] q_1+\boldsymbol{H}[2,2]q_2)+\boldsymbol{H}[3,3]q_3&=&0\\ \boldsymbol{C} [q_1,q_2]^{\mathrm{T}} + \boldsymbol{D} q_3 &=&0\\ q_3&=&-\frac{1}{\boldsymbol{D}} \boldsymbol{C} [q_1,q_2]^{\mathrm{T}} \end{eqnarray}$$

ここで，以下のように表現行列を分解した． $$\boldsymbol{H}=\left( \begin{array}{c} \boldsymbol{A}&\boldsymbol{B}\\\boldsymbol{C}&\boldsymbol{D} \end{array} \right)$$ ただし，$\boldsymbol{A}\in \mathbb{R}^{2\mathrm{x}2}$

この表記からの類推で，さらに高次の場合を考えると，以下が成立する(はず)

$$\begin{eqnarray} \boldsymbol{C} \left( \begin{array}{c} q_1\\q_2 \end{array} \right) + \boldsymbol{D} \left( \begin{array}{c} q_3\\q_4\\ \vdots \end{array} \right) &=&0 \\ \left( \begin{array}{c} q_3\\q_4\\\vdots \end{array} \right)&=& -\boldsymbol{D}^{-1} \boldsymbol{C} \left( \begin{array}{c} q_1 \\ q_2 \end{array} \right) \end{eqnarray}$$

1.の手順の目標である$f(q_3,q_4,\cdots) = g(q_1,q_2)$の形にすることができた(しかもかなり扱いやすい形でラッキー)

2.の手順について，直交条件が使えるように楕円方程式を変形する．

といっても，これは楽勝で，$\boldsymbol{q}_1=[q_1,q_2]^\mathrm{T} , \boldsymbol{q}_2=[q_3,q_4,\cdots]^\mathrm{T}$とおけば以下のようになる．

$$\begin{eqnarray} \left( \begin{array}{c} \boldsymbol{q}_1 \\ \boldsymbol{q}_2 \end{array} \right)^\mathrm{T} \left( \begin{array}{c}\boldsymbol{A}&\boldsymbol{B}\\ \boldsymbol{C}&\boldsymbol{D} \end{array} \right) \left( \begin{array}{c} \boldsymbol{q}_1 \\ \boldsymbol{q}_2 \end{array} \right) &=&\chi^2 \\ \boldsymbol{q}_1^\mathrm{T}(\boldsymbol{Aq}_1+\boldsymbol{Bq}_2)+\boldsymbol{q}_2^\mathrm{T}(\boldsymbol{Cq}_1+\boldsymbol{Dq}_2)&=&\chi^2 \end{eqnarray}$$

3.の手順はもう代入するだけ．

直交条件を$\boldsymbol{q}_1=[q_1,q_2]^\mathrm{T} , \boldsymbol{q}_2=[q_3,q_4,\cdots]^\mathrm{T}$の記法を用いて書き直すと以下のようになる．

$$\begin{eqnarray} \boldsymbol{Cq}_1 + \boldsymbol{Dq}_2 &=&0\\ \boldsymbol{q}_2&=&-\boldsymbol{D}^{-1} \boldsymbol{C} \boldsymbol{q}_1 \end{eqnarray}$$

楕円方程式に代入すれば以下のようになる．

$$\begin{eqnarray} \boldsymbol{q}_1^\mathrm{T}(\boldsymbol{Aq}_1+\boldsymbol{Bq}_2)+\boldsymbol{q}_2^\mathrm{T}(\boldsymbol{Cq}_1+\boldsymbol{Dq}_2)&=&\chi^2\\ \iff\boldsymbol{q}_1^\mathrm{T}(\boldsymbol{Aq}_1+\boldsymbol{Bq}_2)&=&\chi^2\\ \iff \boldsymbol{q}_1^\mathrm{T}(\boldsymbol{Aq}_1-\boldsymbol{B}\boldsymbol{D}^{-1} \boldsymbol{C} \boldsymbol{q}_1) &=& \chi^2\\ \iff \boldsymbol{q}_1^\mathrm{T}( \boldsymbol{A}-\boldsymbol{B}\boldsymbol{D}^{-1} \boldsymbol{C} )\boldsymbol{q}_1 &=& \chi^2 \end{eqnarray}$$

以上より，射影した楕円の表現行列は$\boldsymbol{A}-\boldsymbol{B}\boldsymbol{D}^{-1} \boldsymbol{C}$で与えられることが分かる．